clear all;
close all;
delete(gcp('nocreate'));
%tic
global usual_delta;
usual_delta = 0.05;
% Current folder for Matlab codes
fp.matlab = pwd;
common_path = extractBefore(fp.matlab,"\matlab_dir");
% Paper folder;
fp.paper = sprintf('%s\\paper', common_path);
fp.robust_end_moving = sprintf('%s\\matlab_dir\\robustness\\endogenous_moving', common_path);
fp.robust_two_CES = sprintf('%s\\matlab_dir\\robustness\\counterfactual_two_CES_calibrated', common_path); 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%fixed parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%We estimate the model at lambda=0.95;
fp.lambda = 0.95;

%model primitives;
fp.ages =9:12;
fp.periods=length(fp.ages);
fp.percentiles_kid = linspace(5,95,10);
fp.hc_points_kid = length(fp.percentiles_kid);

fp.Min_Time = 0.01;
fp.Max_Time= 1;
fp.inv_points=20;
fp.percentiles_inv= linspace(5,95,fp.inv_points);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%indexes for latent data
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

fp.ind_data.child_id = 1;
fp.ind_data.child_age = 2;
fp.ind_data.n_sim = 3;
fp.ind_data.child_C = 4;
fp.ind_data.child_C_tp1 = 5;
fp.ind_data.inv = 6;
fp.ind_data.peers = 7;
fp.ind_data.peers_log = 8;
fp.ind_data.ParStyle = 9;
fp.ind_data.peers_tp1 = 10;
fp.ind_data.peers_log_tp1 = 11;
fp.ind_data.neighborhood = 12;
fp.ind_data.moved = 13;

fp.ind_data_network.child_id=1;
fp.ind_data_network.child_age=2;
fp.ind_data_network.n_sim=3;
fp.ind_data_network.child_C=4;
fp.ind_data_network.n_friends=5;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Loading Estimated Moments;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

temp_mean_PS_between_age=importdata('mean_PS_between_class.txt');    
fp.mom_data.mean_PS_between_age = temp_mean_PS_between_age;

moments_temp=importdata('moments.txt');    

variance_moments_temp = importdata('variance_moments.txt');    
fp.variance_moments = variance_moments_temp;

fp.mom_data.mean_PS = moments_temp(1);

fp.ind_mom.coef_PS.child_C = 1;
fp.ind_mom.coef_PS.peers   = 2;
fp.mom_data.coef_PS     = moments_temp(2:3);

fp.ind_mom.coef_skills_tp1.child_C    = 1;
fp.ind_mom.coef_skills_tp1.peers      = 2;
fp.ind_mom.coef_skills_tp1.PS         = 3;

fp.mom_data.coef_skills_tp1 = moments_temp(4:6); %Table 11 Column (1) 
fp.mom_data.coef_skills_tp1_ps0 = moments_temp(7:8);  
fp.mom_data.coef_skills_tp1_ps1 = moments_temp(9:10);

fp.ind_mom.coef_peers_tp1.child_C = 1;
fp.ind_mom.coef_peers_tp1.peers   = 2;
fp.ind_mom.coef_peers_tp1.PS      = 3;

fp.mom_data.coef_peers_tp1  = moments_temp(11:13); %Table 2 Column (2) 
fp.mom_data.coef_peers_tp1_ps0  = moments_temp(14:15); 
fp.mom_data.coef_peers_tp1_ps1  = moments_temp(16:17); 


fp.ind_mom.coef_inv.child_C = 1;
fp.ind_mom.coef_inv.peers   = 2;
fp.ind_mom.coef_inv.const   = 3;

fp.mom_data.coef_inv_PS0  = moments_temp(18:19); 
fp.mom_data.coef_inv_PS1  = moments_temp(20:21); 

fp.mom_data.mean_skills = moments_temp(22:25);

fp.mom_data.mean_inv    = moments_temp(26:27);

fp.mom_data.mean_nfriends = moments_temp(28);


fp.ind_mom.coef_PS_longit.child_C = 1;
fp.ind_mom.coef_PS_longit.peers   = 2;

fp.mom_data.PS_longit = moments_temp(29:30);
fp.mom_data.Inv_longit = moments_temp(31:32);

fp.mom_data.mean_PS_between = moments_temp(33:36);

%Store vector of Moments;
fp.data_mom = moments_temp;

%Loading Initial Conditions (Table 4);
initial_conditions = importdata('initial_conditions.txt');    
fp.schools_distrib = initial_conditions;

%Loading Income information of neighborhoods (it affects preferences for parenting in the model);
mean_income_by_neighborhoods = importdata('family_income_neighborhoods.txt');
mean_income_by_neighborhoods = mean_income_by_neighborhoods./10000;
fp.mean_income_by_neighborhoods = mean_income_by_neighborhoods;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%random mumber draws
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

rng(10); %fix seed;
fp.n_rip_sim = 200;                         % # of simulation of the economy;
fp.tot_draws_mcintegration_network = 50 ;   % # of simulation of the network;

%Initial skills;
fp.init_draws{1} = norminv(rand(fp.schools_distrib(1,3),1,fp.n_rip_sim));
fp.init_draws{2} = norminv(rand(fp.schools_distrib(2,3),1,fp.n_rip_sim));
fp.init_draws{3} = norminv(rand(fp.schools_distrib(3,3),1,fp.n_rip_sim));
fp.init_draws{4} = norminv(rand(fp.schools_distrib(4,3)+50,1,fp.n_rip_sim));

%Peer shocks (backward induction)
fp.peers_backward{1} = rand(length(fp.percentiles_kid)*length(fp.percentiles_kid)*fp.inv_points,fp.schools_distrib(1,3),fp.tot_draws_mcintegration_network, fp.periods-1);
fp.peers_backward{2} = rand(length(fp.percentiles_kid)*length(fp.percentiles_kid)*fp.inv_points,fp.schools_distrib(2,3),fp.tot_draws_mcintegration_network, fp.periods-1);
fp.peers_backward{3} = rand(length(fp.percentiles_kid)*length(fp.percentiles_kid)*fp.inv_points,fp.schools_distrib(3,3),fp.tot_draws_mcintegration_network, fp.periods-1);
fp.peers_backward{4} = rand(length(fp.percentiles_kid)*length(fp.percentiles_kid)*fp.inv_points,fp.schools_distrib(4,3)+50,fp.tot_draws_mcintegration_network, fp.periods-1);

%Peer shocks (forward induction)
fp.peers_forward{1} = rand(fp.schools_distrib(1,3),fp.schools_distrib(1,3),fp.n_rip_sim,fp.periods-1);
fp.peers_forward{2} = rand(fp.schools_distrib(2,3),fp.schools_distrib(2,3),fp.n_rip_sim,fp.periods-1);
fp.peers_forward{3} = rand(fp.schools_distrib(3,3),fp.schools_distrib(3,3),fp.n_rip_sim,fp.periods-1);
fp.peers_forward{4} = rand(fp.schools_distrib(4,3)+50,fp.schools_distrib(4,3)+50,fp.n_rip_sim,fp.periods-1);

%Parenting Style shocks
fp.PS_forward{1} = rand(fp.schools_distrib(1,3),fp.n_rip_sim,fp.periods-1);
fp.PS_forward{2} = rand(fp.schools_distrib(2,3),fp.n_rip_sim,fp.periods-1);
fp.PS_forward{3} = rand(fp.schools_distrib(3,3),fp.n_rip_sim,fp.periods-1);
fp.PS_forward{4} = rand(fp.schools_distrib(4,3)+50,fp.n_rip_sim,fp.periods-1);
    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%COUNTERFACTUALS;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%Set sending and receiving neighborhoods;
fp.origin_neigh         = 1;
fp.receiving_neigh      = 4;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Simulate: Baseline Economy        
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

load('param_est');   
fp.do_estimation =0;
fp.do_counter = 0;            
fp.do_counter_type = 0 ;

obj_eval = obj_function( param_est, fp );
data_baseline_origin = obj_eval.simul_data{fp.origin_neigh};

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;
%Counterfactuals Analysis;

%Figure 3: Treatment Effects of Moving to a Better Neighborhood
%Figure 4: Scaling of Treatment Effects on Skills
%Figure 5: Scaling of Treatment Effects on Parental Behavior
%Figure 6: Endogenous Parental Behavior and Policy Effects

% Table 8: Counterfactual Policy Experiments: Changing Initial Conditions (e.g., Early Childhood Interventions)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;


%Here we double the shocks for peers formation in the backward problem
%because in the counterfactuals we double the state space (moving families
%with different preferences into the new neighborhood);
fp.peers_backward_original = fp.peers_backward;

for j = 1:1:4
   
    fp.peers_backward{j} = [ fp.peers_backward{j} ; fp.peers_backward{j} ] ; 
 
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;
%Policy #1: Treatment Effects of moving a child;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;

fp.do_counter = 1;            
fp.do_counter_type     = 1;
fp.n_kids_moved = 1;

cd(fp.matlab)
%Compute the effect;
obj_counter1 = obj_function( param_est, fp );
%Collect the results;
data_counter1 = obj_counter1.simul_data;
fp.id_moved_children = obj_counter1.id_moved_children;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;
%Figure 3: Treatment Effects of Moving to a Better Neighborhood
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;

policy_effects_counter1(data_counter1,data_baseline_origin,fp)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;
%Policy #2: Policy Effects of moving average many children from Sending Neighborhood);  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;
fp.max_kids_moved  = 90;
fp.min_kids_moved  = 41;
fp.do_counter = 1;            
fp.do_counter_type = 2;

pe_skills_moved      = NaN(1,fp.max_kids_moved);
pe_PS_moved          = NaN(1,fp.max_kids_moved);
pe_inv_moved         = NaN(1,fp.max_kids_moved);
pe_skills_receiving  = NaN(1,fp.max_kids_moved);
pe_PS_receiving      = NaN(1,fp.max_kids_moved);
pe_inv_receiving     = NaN(1,fp.max_kids_moved);

data_baseline = obj_eval.simul_data;


set_children_to_move = NaN(length(fp.min_kids_moved:1:fp.max_kids_moved),fp.n_rip_sim);

for h = 1 : fp.n_rip_sim

set_children_to_move(:,h) = randsample(fp.min_kids_moved:1:fp.max_kids_moved,length(fp.min_kids_moved:1:fp.max_kids_moved));

end


%Loop over the number of children moved. For each policy, we simulate a new
%counterfactual equilibrium of the economy;
ind = 1;
for n_kids =  [1, 5:5:50]
    
    n_kids
    
    fp.n_kids_moved    = n_kids;

    fp.shock_child_moved = set_children_to_move(1:n_kids,:);

    cd(fp.matlab)
    obj_counter2 = obj_function( param_est, fp );
    
    data_counter2 = obj_counter2.simul_data;
    fp.id_moved_children = obj_counter2.id_moved_children;

    pe = policy_effects_counter2(data_counter2,data_baseline,fp);
    
    pe_skills_moved(ind)=pe.pe_skills_moved;
	pe_PS_moved(ind) = pe.pe_PS_moved;
	pe_inv_moved(ind) = pe.pe_inv_moved;
    pe_skills_receiving(ind) = pe.pe_skills_receiving;
	pe_PS_receiving(ind) = pe.pe_PS_receiving;
    pe_inv_receiving(ind) = pe.pe_inv_receiving;
           
    ind = ind + 1;
end



grid_x= [1, 5:5:50];
grid_y = grid_x;
cd(fp.paper)


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;
%Figure 4: Scaling of Treatment Effects on Skills
%Figure 5: Scaling of Treatment Effects on Parental Behavior
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;

h1=figure;
scatter(grid_x ,pe_skills_moved(1:length(grid_y)) )
hold on
p = polyfit(grid_x ,pe_skills_moved(1:length(grid_y)),2);
y = polyval(p, grid_y );
plot( grid_x , y , '--r','LineWidth',2,'MarkerSize',10)
legend('Policy Effects','Quadratic Approx')
xlabel('N. Children Moved');
ylabel('Mean Effect on Log-Skills (Moved Children)');
set(gca,'XTick',0:5:length(fp.min_kids_moved:1:fp.max_kids_moved))
set(gca,'YTick',0.06:0.02:0.24)
set(gca,'YLim', [0.06 0.24 ])
%hgexport(h1, 'Count1.png',hgexport('factorystyle'), 'Format', 'png');
%hgexport(h1, 'Count1.eps',hgexport('factorystyle'), 'Format', 'eps');
    

h2=figure;
scatter(grid_x ,pe_PS_moved(1:length(grid_y)) )
hold on
p = polyfit(grid_x ,pe_PS_moved(1:length(grid_y)),2);
y = polyval(p,grid_y );
plot( grid_x , y , '--r','LineWidth',2,'MarkerSize',10)  
plot( grid_x , zeros(length(grid_x),1),'--k','LineWidth',1,'MarkerSize',3)
legend('Policy Effects','Quadratic Approx','Location','NorthEast')
xlabel('N. Children Moved');
ylabel('Mean Effect on Author. Parenting (Moved Children)');
axis([0 length(fp.min_kids_moved:1:fp.max_kids_moved) -0.2 0.05])
%hgexport(h2, 'Count2.png',hgexport('factorystyle'), 'Format', 'png'); 
%hgexport(h2, 'Count2.eps',hgexport('factorystyle'), 'Format', 'eps'); 

h3=figure;
scatter( grid_x ,pe_inv_moved(1:length(grid_y)) )
hold on
p = polyfit( grid_x ,pe_inv_moved(1:length(grid_y)),2);
y = polyval(p, grid_y );
plot( grid_x , y , '--r','LineWidth',2,'MarkerSize',10)
plot( grid_x , zeros(length(grid_x),1),'--k','LineWidth',1,'MarkerSize',3)    
legend('Policy Effects','Quadratic Approx')    
xlabel('N. Children Moved');
ylabel('Mean Effect on Time Investment (Moved Children)');
axis([0 length(fp.min_kids_moved:1:fp.max_kids_moved) -0.1 0.055])
%hgexport(h3, 'Count3.png',hgexport('factorystyle'), 'Format', 'png');
%hgexport(h3, 'Count3.eps',hgexport('factorystyle'), 'Format', 'eps');


h4=figure;
scatter( grid_x ,pe_skills_receiving(1:length(grid_y)) )
hold on
p = polyfit( grid_x ,pe_skills_receiving(1:length(grid_y)),2);
y = polyval(p, grid_y );
plot( grid_x , y , '--r','LineWidth',2,'MarkerSize',10)
legend('Policy Effects','Quadratic Approx')
xlabel('N. Children Moved');
ylabel('Mean Effect on Log-Skills (Receiving Children)');
%hgexport(h4, 'Count4.png',hgexport('factorystyle'), 'Format', 'png');    
%hgexport(h4, 'Count4.eps',hgexport('factorystyle'), 'Format', 'eps');    

h5=figure;
scatter( grid_x ,pe_PS_receiving(1:length(grid_y)) )
hold on
p = polyfit( grid_x ,pe_PS_receiving(1:length(grid_y)),2);
y = polyval(p, grid_y );
plot(grid_x , y , '--r','LineWidth',2,'MarkerSize',10)
plot( grid_x , zeros(length(grid_x),1),'--k','LineWidth',1,'MarkerSize',3)    
legend('Policy Effects','Quadratic Approx','Location','NorthWest')      
xlabel('N. Children Moved');
ylabel('Mean Effect on Author. Parenting (Receiving Children)');
axis([0 length(fp.min_kids_moved:1:fp.max_kids_moved) -0.01 0.02])
%hgexport(h5, 'Count5.png',hgexport('factorystyle'), 'Format', 'png'); 
%hgexport(h5, 'Count5.eps',hgexport('factorystyle'), 'Format', 'eps'); 

h6=figure;
scatter( grid_x ,pe_inv_receiving(1:length(grid_y)) )
hold on
p = polyfit( grid_x ,pe_inv_receiving(1:length(grid_y)),2);
y = polyval(p, grid_y );
plot(grid_x , y , '--r','LineWidth',2,'MarkerSize',10)
plot( grid_x , zeros(length(grid_x),1),'--k','LineWidth',1,'MarkerSize',3)    
legend('Policy Effects','Quadratic Approx','Location','NorthWest') 
xlabel('N. Children Moved');
ylabel('Mean Effect on Time Investment (Receiving Children)');
axis([0 length(fp.min_kids_moved:1:fp.max_kids_moved) -0.01 0.05])
%hgexport(h6, 'Count6.png',hgexport('factorystyle'), 'Format', 'png');   
%hgexport(h6, 'Count6.eps',hgexport('factorystyle'), 'Format', 'eps');   
    
cd(fp.matlab)


%Store policy results for skills;
save('pe_skills_moved','pe_skills_moved')
save('pe_skills_receiving','pe_skills_receiving')

%Save in this folder to compare these results with alternative models;
cd(fp.robust_end_moving)
save('pe_skills_moved','pe_skills_moved')
save('pe_skills_receiving','pe_skills_receiving')

cd(fp.robust_two_CES)
save('pe_skills_moved','pe_skills_moved')
save('pe_skills_receiving','pe_skills_receiving')
cd(fp.matlab)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;
%Policy #3: Policy #2 with fixed Parental Behavior;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;
fp.do_counter          = 1;
fp.do_counter_type     = 3;

pe_skills_moved_fixparenting      = NaN(1,fp.max_kids_moved);
pe_skills_receiving_fixparenting  = NaN(1,fp.max_kids_moved);

fp.shock_child_moved = [];
fp.simul_data_baseline = data_baseline;



%Loop over the number of children moved, keeping parental behavior fixed at
%the baseline level;
ind = 1;
for n_kids =  [1, 5:5:50]
    
    n_kids
    
    fp.n_kids_moved    = n_kids;

    fp.shock_child_moved = set_children_to_move(1:n_kids,:);
    
    cd(fp.matlab)
    obj_eval_fixparenting = obj_function( param_est, fp );
    
    simul_data_counterfactual_fixparenting = obj_eval_fixparenting.simul_data;    
    fp.id_moved_children = obj_eval_fixparenting.id_moved_children;
        
    pe = policy_effects_counter2(simul_data_counterfactual_fixparenting,data_baseline,fp);
    
    pe_skills_moved_fixparenting(ind)=pe.pe_skills_moved;
    pe_skills_receiving_fixparenting(ind) = pe.pe_skills_receiving;
    
    ind = ind + 1 ;
    
    
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;
%Figure 6: Endogenous Parental Behavior and Policy Effects
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;

    cd(fp.paper)

    h8=figure;
    scatter(grid_x ,pe_skills_moved(1:length(grid_y)) ,  'b')
    hold on
    p = polyfit(grid_x ,pe_skills_moved(1:length(grid_y)),2);
    y = polyval(p, grid_y );
    plot( grid_x , y , '--b','LineWidth',2,'MarkerSize',10)
    
    scatter(grid_x ,pe_skills_moved_fixparenting(1:length(grid_y)) ,  'r')
    p = polyfit(grid_x ,pe_skills_moved_fixparenting(1:length(grid_y)),2);
    y = polyval(p, grid_y );
    plot( grid_x , y , '--r','LineWidth',2,'MarkerSize',10)    
    
    legend('Policy Effects','Quadratic Approx', 'Policy Effects (Fixed Parent Behavior)' , 'Quadratic Approx' )
    xlabel('N. Children Moved');
    ylabel('Mean Effect on Log-Skills (Moved Children)');
    set(gca,'XTick',0:5:length(fp.min_kids_moved:1:fp.max_kids_moved))
    set(gca,'YTick',0.08:0.02:0.28)
    set(gca,'YLim', [0.08 0.28])
    %hgexport(h8, 'Count_fix1.png',hgexport('factorystyle'), 'Format', 'png');     
    %hgexport(h8, 'Count_fix1.eps',hgexport('factorystyle'), 'Format', 'eps');     
    

    h9=figure;
    scatter( grid_x ,pe_skills_receiving(1:length(grid_y)) ,  'b')
    hold on
    p = polyfit( grid_x ,pe_skills_receiving(1:length(grid_y)),2);
    y = polyval(p, grid_y );
    plot( grid_x , y , '--b','LineWidth',2,'MarkerSize',10)
    
    scatter( grid_x ,pe_skills_receiving_fixparenting(1:length(grid_y)) ,  'r')
    p = polyfit( grid_x ,pe_skills_receiving_fixparenting(1:length(grid_y)),2);
    y = polyval(p, grid_y );
    plot( grid_x , y , '--r','LineWidth',2,'MarkerSize',10)             
    legend('Policy Effects','Quadratic Approx', 'Policy Effects (Fixed Parent Behavior)' , 'Quadratic Approx' )
    xlabel('N. Children Moved');
    ylabel('Mean Effect on Log-Skills (Receiving Children)');
    set(gca,'YTick',-0.18:0.02:0.00)
    set(gca,'YLim', [-0.18 0])    
    %hgexport(h9, 'Count_fix2.png',hgexport('factorystyle'), 'Format', 'png');
    %hgexport(h9, 'Count_fix2.eps',hgexport('factorystyle'), 'Format', 'eps');           
    


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Other Set of Counterfactual Analysis/Decompositions;

% Table 8: Counterfactual Policy Experiments: Changing Initial Conditions (e.g., Early Childhood Interventions)

% 1. No Inequality(all neighborhoods share initial aggregate mean)
% 2. No Between-Neighb. Inequality (all initial conditions = aggreg mean and SD)
% 3. No Within-Neighb. Inequality 
% 4. Truncate Local Distrib. at 10th percent (truncating within each neighborhod)
% 5. Halving Cost of Parental Investments
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
cd(fp.matlab)
fp.do_counter = 1;            

fp.peers_backward = fp.peers_backward_original;

data_baseline = obj_eval.simul_data;

data_baseline_stacked = [];


for j = 1:4
    
    data_baseline_stacked = [data_baseline_stacked ;  log(data_baseline{j}(  data_baseline{j}(:,fp.ind_data.child_age , 1 )==9, fp.ind_data.child_C , 1 ))    ];
    
end

fp.national_mean = nanmean(data_baseline_stacked);
fp.national_sd   = nanstd(data_baseline_stacked);


%%%%%%%%%%%%%%%%%%%%%%%%%%%;
% 1. all location have the national mean, but no within-inequality (No Inequality);
%%%%%%%%%%%%%%%%%%%%%%%%%%%;    
    
fp.do_counter_type = 4;    
obj_no_inequality = obj_function( param_est, fp );
data_no_inequality = obj_no_inequality.simul_data;

%%%%%%%%%%%%%%%%%%%%%%%%%%%;
% 2. all location have the same initial mean and SD, equal to the national average (No Between-Neighborhood Inequality);
%%%%%%%%%%%%%%%%%%%%%%%%%%%;     

fp.do_counter_type = 5;    
obj_no_between_inequality = obj_function( param_est, fp );
data_no_between_inequality = obj_no_between_inequality.simul_data;

%%%%%%%%%%%%%%%%%%%%%%%%%%%;
% 3. all location have the original mean, but no within-inequality (No Within-Neighborhood Inequality) ;
%%%%%%%%%%%%%%%%%%%%%%%%%%%;

fp.do_counter_type = 6;    
obj_no_within_inequality = obj_function( param_est, fp );
data_no_within_inequality = obj_no_within_inequality.simul_data;

%%%%%%%%%%%%%%%%%%%%%%%%%%%;
% 4. Truncate Local Distrib. at 10th percent (truncating within each neighborhod)
%%%%%%%%%%%%%%%%%%%%%%%%%%%;

fp.do_counter_type = 7;
obj_truncated_10th = obj_function( param_est, fp );
data_truncated_10th = obj_truncated_10th.simul_data;

%%%%%%%%%%%%%%%%%%%%%%%%%%%;
% 5. Halving Cost of Parental Investments
%%%%%%%%%%%%%%%%%%%%%%%%%%%;

fp.do_counter_type = 8;
obj_halved_cost_inv = obj_function( param_est, fp );
data_halved_cost_inv = obj_halved_cost_inv.simul_data;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;
% Table 8: Counterfactual Policy Experiments: Changing Initial Conditions (e.g., Early Childhood Interventions)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;

New_Counter_Tables(data_baseline, data_no_inequality, data_no_between_inequality, data_no_within_inequality, data_truncated_10th, data_halved_cost_inv , fp)    


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;
% How Policy Effects Vary with Structural Parameters;
% %Policy # 2 with 50 Children moved. For each new parameter value:
    %1. We compute the new baseline economy;
    %2. We compute the new counterfactual economy (with 50 children moved)
    %3. We compare baseline vs. counterfactual;

%Appendix Figure E-2: The Impact of Homophily on the Policy Scale Up
%Appendix Figure E-3: Parental Responses and the Elasticity of Substitution Peers vs. Parents

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;

%Technology (Cobb-Douglas, PS=1);
param_con(1) = exp(param_est(1)/100);    %Children's skills (alpha1);
param_con(2) = exp(param_est(2)/100);    %Peers' Skills (alpha2);
param_con(3) = exp(param_est(3)/100);    %Investments (alpha3);

param_con(4) = param_est(4)/100;         %Parenting Style (alpha0);

%Technology (CES, PS=0);
param_con(5) = exp(param_est(5)/100)./(1 + exp(param_est(5)/100) );         %CES Elasticity (   Parents vs Peers, inner)
param_con(6) = exp(param_est(6)/100)./(1 + exp(param_est(6)/100) );         %CES Share Skills ;
param_con(7) = exp(param_est(7)/100)./(1 + exp(param_est(7)/100) );         %CES Share Peers;
param_con(8) = -exp(param_est(8)/100) ;                                     %CES Elasticity ( Skills vs Parents-Peers, outer);
param_con(9) = param_est(9)/100 ;                                           %Return to Scale CES;

param_con(10) = param_est(10)/100;           %TFP constant;
param_con(11) = param_est(11)/100;           %TFP trend;

ind  = 11 ;

%Parent's Preferences;

param_con(ind+1) =      exp(param_est(ind+1)/100);                    %Weight on Children's Skills;
param_con(ind+2) =      param_est(ind+2)/100 ;                        %Disutility of PS;
param_con(ind+3) =      param_est(ind+3)/100 ;                        %Heterogeneity in PS;

ind  = ind + 3 ;

%Child's Preferences;

param_con(ind+1) =  param_est(ind+1)/100;                       %Constant (gamma0);
param_con(ind+2) =  param_est(ind+2)/100;                       %Own Skills (gamma1);
param_con(ind+3) =  param_est(ind+3)/100;                       %Child j 's Skills (gamma2);
param_con(ind+4) =  param_est(ind+4)/100;                       %Homophily (gamma3);
param_con(ind+5) =  param_est(ind+5)/100;                       %PS effect (gamma4);

param_con(ind+6) =      param_est(ind+6)/100 ;                  %Disutility of PS (by Neighborhood Income);


grid_perturb = 0.5:0.1:1.5;

pe_skills_moved_perturb1     = NaN(1,length(grid_perturb));
pe_skills_receiving_perturb1  = NaN(1,length(grid_perturb));

param_con_pertub = param_con;

% We move 50 children;    
ind = 50;
fp.n_kids_moved    = ind;
fp.shock_child_moved = set_children_to_move(1:ind,:);



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;
%Appendix Figure E-2: The Impact of Homophily on the Policy Scale Up
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;

for j = 1:1:length(grid_perturb)
param_con_pertub(18) = param_con(18).*grid_perturb(j);
cd(fp.matlab)

%%%%%%%%%%%%%%%%%%%;
%Baseline;
%%%%%%%%%%%%%%%%%%%;

fp.do_counter = 0;            
fp.do_counter_type = 0 ;

%reset shocks;
fp.peers_backward = fp.peers_backward_original;
obj_base_perturb = sim_data_function( param_con_pertub, fp );
data_baseline_perturb = obj_base_perturb.sim_data;


%%%%%%%%%%%%%%%%%%%;
%Counterfactual;
%%%%%%%%%%%%%%%%%%%;

fp.do_counter = 1;            
fp.do_counter_type = 2;

%Here we double the shocks for peers formation in the backward problem
%because in the counterfactuals we double the state space (moving families
%with different preferences into the new neighborhood);
for i = 1:1:4
   
    fp.peers_backward{i} = [ fp.peers_backward_original{i} ; fp.peers_backward_original{i} ] ; 
 
end

obj_counter2_perturb = sim_data_function(param_con_pertub,fp) ;
data_counter2_perturb = obj_counter2_perturb.sim_data ; 
fp.id_moved_children = obj_counter2_perturb.id_moved_children;

pe = policy_effects_counter2(data_counter2_perturb,data_baseline_perturb,fp);

pe_skills_moved_perturb1(j)=pe.pe_skills_moved;

pe_skills_receiving_perturb1(j) = pe.pe_skills_receiving;
end



h10=figure;
plot(grid_perturb,pe_skills_moved_perturb1(1,:),'--bo','LineWidth',3.5,...
                       'MarkerEdgeColor','b',...
                       'MarkerFaceColor','b',...
                       'MarkerSize',7.5)
hold on
xlabel('Homophily Parameter (\gamma_3, relative to baseline estimates)');
ylabel('Policy Effects on Moved Children');
hold off
cd(fp.paper)    
%hgexport(h10, 'Policy_Perturbation1.png',hgexport('factorystyle'), 'Format', 'png');
%hgexport(h10, 'Policy_Perturbation1.eps',hgexport('factorystyle'), 'Format', 'eps'); 
cd(fp.matlab)



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;
%Appendix Figure E-3: Parental Responses and the Elasticity of Substitution Peers vs. Parents
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;

grid_perturb_alpha30 = grid_perturb(3:end-2);

%reset initial baseline estimates;
param_con_pertub = param_con;

pe_skills_moved_perturb2    = NaN(1,length(grid_perturb_alpha30));
pe_skills_receiving_perturb2  = NaN(1,length(grid_perturb_alpha30));

pe_skills_receiving_perturb2_base = NaN(1,length(grid_perturb_alpha30));
pe_skills_receiving_perturb2_count = NaN(1,length(grid_perturb_alpha30));

pe_inv_moved_perturb2  = NaN(1,length(grid_perturb_alpha30));
pe_inv_receiving_perturb2  = NaN(1,length(grid_perturb_alpha30));

for j = 1:1:length(grid_perturb_alpha30)
param_con_pertub(5) = param_con(5).*grid_perturb_alpha30(j);

%To prevent that the CES elasticity parameter goes above 1;
if param_con_pertub(5)>1
param_con_pertub(5) = 1;
end

%%%%%%%%%%%%%%%%%%%;
%Baseline;
%%%%%%%%%%%%%%%%%%%;

fp.do_counter = 0;            
fp.do_counter_type = 0 ;

%reset shocks;
fp.peers_backward = fp.peers_backward_original;
obj_base_perturb = sim_data_function( param_con_pertub, fp );
data_baseline_perturb = obj_base_perturb.sim_data;

%%%%%%%%%%%%%%%%%%%;
%Counterfactual;
%%%%%%%%%%%%%%%%%%%;


fp.do_counter = 1;            
fp.do_counter_type = 2;

%Here we double the shocks for peers formation in the backward problem
%because in the counterfactuals we double the state space (moving families
%with different preferences into the new neighborhood);
for i = 1:1:4
   
    fp.peers_backward{i} = [ fp.peers_backward_original{i} ; fp.peers_backward_original{i} ] ; 
 
end

cd(fp.matlab)
obj_counter2_perturb = sim_data_function(param_con_pertub,fp) ;
data_counter2_perturb = obj_counter2_perturb.sim_data ; 
fp.id_moved_children = obj_counter2_perturb.id_moved_children;

pe = policy_effects_counter2(data_counter2_perturb,data_baseline_perturb,fp);


pe_inv_moved_perturb2(j)     = pe.pe_inv_moved;  
pe_inv_receiving_perturb2(j) = pe.pe_inv_receiving;
end




h12=figure;
plot(grid_perturb_alpha30,pe_inv_receiving_perturb2,'--bo','LineWidth',3.5,...
                       'MarkerEdgeColor','b',...
                       'MarkerFaceColor','b',...
                       'MarkerSize',7.5)
xlabel('Elast. Substitution Peers vs. Parents (\alpha_{3,0}, relative to baseline estimates)');
ylabel('Policy Effects on Investments of Receiving Families');
xlim([0.7 1.3])
cd(fp.paper)    
%hgexport(h12, 'Policy_Perturbation2_inv_receiving.png',hgexport('factorystyle'), 'Format', 'png');
%hgexport(h12, 'Policy_Perturbation2_inv_receiving.eps',hgexport('factorystyle'), 'Format', 'eps'); 
cd(fp.matlab)


h13=figure;
plot(grid_perturb_alpha30,pe_inv_moved_perturb2,'--bo','LineWidth',3.5,...
'MarkerEdgeColor','b',...
'MarkerFaceColor','b',...
'MarkerSize',7.5)
xlabel('Elast. Substitution Peers vs. Parents (\alpha_{3,0}, relative to baseline estimates)');
ylabel('Policy Effects on Investments of Moved Families');
xlim([0.7 1.3])
cd(fp.paper)
%hgexport(h13, 'Policy_Perturbation2_inv_moved.png',hgexport('factorystyle'), 'Format', 'png');
%hgexport(h13, 'Policy_Perturbation2_inv_moved.eps',hgexport('factorystyle'), 'Format', 'eps'); 
cd(fp.matlab)
